clear all;
%% Code explanation
% 170305 DZ 
% 

%% Define the resistance to temperature conversion equation
T_RuOx=@(R)1./exp(-1.001200439904+1.175279377873*log(R-2.2)+0.021010648114*(log(R-2.2)).^2-0.024007137972*(log(R-2.2)).^3+0.002054398372*(log(R-2.2)).^4);
%
%% Set the parameter

L_W1=121.9/105.5; % 3.3Sn/12PbTe #766
L_W2=85.2/105.5; % 3Sn/12PbTe #849

ratio=0.5;
%% Read the data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

root=('D:\Documents2\Rawdata\Sn\171219 MPI mK\20180118 3Sn_12PbTe out_of_plane_data\'); %

filenames={'20180118_3samples_1200uW_v6145p_onlyrotary_perp1';'20180118_3samples_850uW_v6145p_onlyrotary_perp1';...
    '20180118_3samples_600uW_v6145p_onlyrotary_perp1';'20180118_3samples_500uW_v615p_onlyrotary_perp1';'20180118_3samples_400uW_v6155p_onlyrotary_perp1';...
    '20180118_3samples_300uW_v6155p_onlyrotary_perp1';'20180118_3samples_200uW_v6155p_onlyrotary_perp1';'20180118_3samples_200uW_v6155p_onlyrotary_perp2';...
    '20180118_3samples_0uW_v6155p_onlyrotary_perp1';'20180119_3samples_0uW_v6155p_onlyrotary_perp2';'20180119_3samples_0uW_v6162p_onlyrotary_perp1';...
    '20180119_3samples_0uW_v6175p_onlyrotary_perp1';'20180119_3samples_0uW_v622p_onlyrotary_perp1';'20180119_3samples_0uW_v630p_onlyrotary_perp1';...
    '20180119_3samples_0uW_v630p_onlyrotary_perp2';'20180119_3samples_0uW_v630p_onlyrotary_perp3';'20180119_3samples_0uW_v6100p_onlyrotary_perp3';...
    '20180120_3samples_250uW_v6100p_withroots_perp1';'20180120_3samples_180uW_v6100p_withroots_perp1';'20180120_3samples_80uW_v6100p_withroots_perp1';...
    '20180120_3samples_50uW_v6100p_withroots_perp1';'20180120_3samples_20uW_v6100p_withroots_perp1';'20180121_3samples_0uW_v6100p_withroots_perp1'};



Boffset=0.009;
NT=length(filenames);


colorset=jet(NT);
figure

Tmc=zeros(NT,1);
Tmc_err=zeros(NT,1);
Hc2=zeros(NT,2);

kk=NT-3;
file1=[root,filenames{kk},'_measdown.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%s %s %s %s %s %s','commentStyle','#');
        data1=data1';
        for ii=1:size(data1,1)
            data1{ii,:}=str2double(data1{ii,:});
        end       
        fclose(fid);
   
          B=data1{1}-Boffset;
          R1=(data1{2}/L_W1+0.009)/5*1000; % factor 5 comes because the current was 50 nA, lock-in sens: 100 uV
          R2=(data1{3}/L_W2+0.02)/5*1000;
indexB=find(B>1.4);
p1=polyfit(B(indexB),R1(indexB),1);
p2=polyfit(B(indexB),R2(indexB),1);
 subplot(1,2,1),plot(B,ratio*(p1(1)*B+p1(2)),'-k');
          hold on
 subplot(1,2,2),plot(B,ratio*(p2(1)*B+p2(2)),'-k');
          hold on

for kk=1:NT
    jj=NT-kk+1;
    file1=[root,filenames{jj},'_measdown.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%s %s %s %s %s %s','commentStyle','#');
        data1=data1';
        for ii=1:size(data1,1)
            data1{ii,:}=str2double(data1{ii,:});
        end       
        fclose(fid);
   
          B=data1{1}-Boffset;
          R1=(data1{2}/L_W1+0.009)/5*1000;
          R2=(data1{3}/L_W2+0.02)/5*1000; 
%           R3=data1{4};
          T=T_RuOx(data1{5}*100);
          Tmc(kk)=mean(T);
          Tmc_err(kk)=3*std(T);
          
          subplot(1,2,1),plot(B,R1,'Color',colorset(kk,:)','LineWidth',1);
          hold on
           subplot(1,2,2),plot(B,R2,'Color',colorset(kk,:)','LineWidth',1);
          hold on

        index1=find((R1-ratio*p1(1)*B-ratio*p1(2)<=0));
        if isempty(index1)==0
          Hc2(kk,1)=B(index1(1));
        else
            Hc2(kk,1)=NaN;
        end
          
          index2=find((R2-ratio*p2(1)*B-ratio*p2(2)<=0));
        if isempty(index2)==0
          Hc2(kk,2)=B(index2(1));
        else
            Hc2(kk,2)=NaN;
        end
        
        title={'B(T)','R(Ohm)','T(K)','averaged T(K)'};
        AA=[B R2 T Tmc(kk)*ones(length(B),1)];
       
         xlswrite('Fig2D.xlsx',AA,kk);
       xlswrite('Fig2D.xlsx',title,kk);
%     
end

index1=find(isnan(Hc2(:,1))==0);
index2=find(isnan(Hc2(:,2))==0);

Hc21=Hc2(index1,1);
Hc22=Hc2(index2,2);
Tmc1=Tmc(index1);
Tmc2=Tmc(index2);
N1=length(Tmc1);
M1=N1-2;
N2=length(Tmc2);
M2=N2-3;

pp1=polyfit(Tmc1(M1:N1),Hc21(M1:N1),1);
pp2=polyfit(Tmc2(M2:N2),Hc22(M2:N2),1);

Tmc10=abs(pp1(2)/pp1(1));
Tmc20=abs(pp2(2)/pp2(1));

xi2=sqrt(6.626E-34/4/pi/1.6E-19/(abs(pp2(2))))*1E9;
Tp=(0:0.002:Tmc20)';
figure
plot([Tmc1' Tmc10]',[Hc21' 0]','dr',[Tmc2' Tmc20]',[Hc22' 0]','sb',Tp,pp2(1)*Tp+pp2(2),'-k');
hold on





%% End




